Lyapunov Modes and Time- Correlation 
Functions for Two-Dimensional Systems 



Tooru Taniguchi and Gary P. Morriss 

School of Physics, University of New South Wales, Sydney, New South Wales 

2052, Australia 



Abstract 

The relation between the Lyapunov modes (delocalized Lyapunov vectors) and 
the momentum autocorrelation function is discussed in two-dimensional hard-disk 
systems. We show numerical evidence that the smallest time-oscillating period of 
the Lyapunov modes is twice as long as the time-oscillating period of momentum 
autocorrelation function for both square and rectangular two-dimensional systems 
with hard-wall boundary conditions. 
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1 Introduction 



The Lyapunov mode is a delocalized structure appearing in the Lyapunov 
vectors corresponding to the Lyapunov exponents close to zero in many-body 
chaotic systems. The existence of such a mode structure for Lyapunov vec- 
tors was first suggested from the stepwise structure of the Lyapunov spectra 
for many-hard-core-particle systems [1,2], and led to the discovery of a sta- 
tionary transverse T-mode structure in the Lyapunov vectors [2,3,4]. Further 
investigations clarified the other two kinds of Lyapunov modes, the longitudi- 
nal L-mode [5] and the momentum proportional P-mode [6,7,8], both of which 
contain an explicit time-dependence. Combining these three kinds of Lyapunov 
modes, explained the degeneracies in the stepwise structure of Lyapunov spec- 
tra [7,8,9]. Moreover, the structure of these Lyapunov modes suggests that the 
origin of the Lyapunov modes is the dynamical conservation laws, like energy 
and momentum conservation, and spatial and temporal translational invari- 
ances [6,7,8,9] and Noether's theorem. The conservation laws and translational 
invariances dominate the global and thermodynamic behavior of many-body 
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systems, and the behavior of Lyapunov vectors for the Lyapunov exponents 
close to zero should correspond to slow and collective movements. Therefore 
it is expected that the Lyapunov modes (which are chaotic properties) allow 
us to discuss the thermodynamic behavior of many-body systems dynami- 
cally. Many studies of the Lyapunov modes have used hard-core-particles in 
quasi-one-, two- or three-dimensional systems because of its established Lya- 
punov vector dynamics and clear mode structures. Recently Lyapunov modes 
for soft-core particle systems have also been investigated [10,11,12]. Some an- 
alytical approaches have been proposed to understand the Lyapunov modes 
and the stepwise structure of Lyapunov spectra. These include a random ma- 
trix approach [13], a Fokker-Planck equation approach [14], kinetic approaches 
[15,16,17] and a periodic orbit approach [18]. 

One recent development in the study of Lyapunov modes is the discovery 
of a connection between the time dependence of the Lyapunov modes and 
the momentum autocorrelation function [7,8]. The integral of the momentum 
autocorrelation function gives the transport coefficient (here the diffusion co- 
efficient) through linear response theory [19]. Moreover it includes information 
on the collective movement of many-body systems [20,21,22], and is accessi- 
ble experimentally using neutron and light scattering techniques [23,24,25]. It 
was shown in Ref. [7,8] that the largest time-oscillating period T^ ya of time- 
dependent Lyapunov modes is twice as long as the time-oscillating period T tc f 
of the momentum autocorrelation function, namely 



Lya 



= 2T te/ . (1) 



This relation connects the Lyapunov mode and an experimentally accessible 
quantity. In reference [7,8], a theoretical argument to support this relation, 
and numerical evidence that this relation is independent of density, particle 
number, and boundary conditions in quasi-one- dimensional systems was given. 
The aim of this paper is to present numerical evidence that the relation (1) 
is satisfied not only in quasi-one-dimensional systems but also in fully two- 
dimensional systems. 

To calculate Lyapunov exponents and Lyapunov vectors we use the Lyapunov 
vector dynamics for hard-disk systems [1], and the algorithm developed by 
Benettin et al. [26], and Shimada and Nagashima [27]. In general, the dif- 
ference between two nearby phase space trajectories defines a direction in 
phase space. The rate of growth of this difference is controlled by the largest 
Lyapunov exponent, and the normalized direction (which is time dependent) 
is the Lyapunov vector for that exponent. Eliminating the direction defined 
by the Lyapunov vector for the largest exponent, reduces the dimension of 
the phase space by one. In that reduced phase space we can again define a 
direction that corresponds to the (next) largest Lyapunov exponent and its 
associated Lyapunov vector. Continuing this procedure we define the same 
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Fig. 1. The stepwise structure of Lyapunov spectrum for a square system consisting 
of 100 hard-disks with hard- wall boundary conditions (circles). For a comparison, 
the Lyapunov spectrum for a square system with periodic boundary conditions 
(crosses) is also shown. The Lyapunov exponent indicated by an arrow is the one in 
which a time-oscillation of Lyapunov mode is shown in Fig. 2(b). Inset: Full positive 
branch of Lyapunov spectra. 



number of Lyapunov exponents and Lyapunov vectors as there are indepen- 
dent directions in phase space. This process is essentially the Benettin scheme 
[26] for calculating the full spectrum of Lyapunov exponents for the system. 
The only additional element is the significance given to basis vectors used in 
the Benettin scheme; here they called Lyapunov vectors and become an object 
of study in their own right. While the Lyapunov vectors associated with the 
largest exponents are strongly localized (having non-zero components for only 
a few particles) the vectors associated with the smallest positive and negative 
exponents have delocalized Lyapunov vectors that we refer to as Lyapunov 
modes. 



2 A Square System with 100 Disks 

We investigate the time-oscillating periods of the Lyapunov modes and the 
momentum autocorrelation function, for a system of hard-disks in a square 
with the hard-wall boundary conditions. Figure 1 is the Lyapunov spectrum 
for this system (circles). Here, we used as system parameters: N = 100, disk 
radius R = 1, mass M = 1, total energy E = N, and the system lengths in 
the x and y directions L x = L y ps 38.7, so that the particle number density is 
N(2R) 2 /(L x L y ) ps 0.267. The Lyapunov spectrum shows a 4-point step close 
to the two zero-Lyapunov exponents. We use hard-wall boundary conditions 
because that leads to more and clearer steps in the Lyapunov spectrum than 
for the equivalent system with periodic boundary conditions. In fact we do not 
observe stepwise structure of the Lyapunov spectrum for a 100-disk system 
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Fig. 2. (a) The spatial region of the system in which the average of the longitudinal 
Lyapunov mode of Fig. (b), is taken, (b) The time-oscillation of the local spatial re- 
gion for the Lyapunov vector component {5y^ 2N ~ 1 ^)t for the component dy^ of 
the Lyapunov vector. The line is a fit of the numerical data to a sinusoidal function. 
The system is the same as that in Fig. 1 with hard-wall boundary conditions. 



with periodic boundary conditions. To show this, in Fig. 1, we also plotted 
the Lyapunov spectrum for a square system with periodic boundary conditions 
(crosses). (Here we used the same system parameters as for hard- wall bound- 
ary conditions, except that L x = L y « 36.7 so the effective region available 
for the disks to move is the same in both cases.) The system with periodic 
boundary conditions is too small to observe the stepwise structure of the Lya- 
punov spectrum. The number of zero-Lyapunov exponents is 2 for hard-wall 
boundary conditions, because of energy conservation and time-translational 
invariance (determinisity of orbit), and is 6 for periodic boundary conditions, 
because of total momentum conservation, conservation of center of mass in 
both directions, energy conservation and time translational invariance [28]. 
Only half of those zero-Lyapunov exponents are presented in Fig. 1. 

We are interested in the time-oscillation of modes corresponding to the step- 
wise structure of the Lyapunov spectrum shown in Fig. 1. It is known that 
there are two kinds of time-dependent Lyapunov modes, the L-modes and the 
P-modes, each with the same period [6,7,8,9]. Our concern is to determine 
the time-oscillating period of Lyapunov mode, so it is enough to consider one 
of these. We choose the L-mode and concentrate into its time-oscillating be- 
havior only, ignoring the spatial structure [7,8,9]. To determine the period of 
the L-mode we use the same system with hard-wall boundary conditions, see 
Fig. 1, and consider the fiyf N . This Lyapunov vector component Sy^ 
corresponds to the Lyapunov exponent A^" 1 -* and is indicated by the arrow 
in Fig. 1. To get clear time-oscillating behavior we take the average of 5yf N ~ lS> 
for disks whose centers (xj,yj) are in the region R < Xj < R + [L x — 2R)/3 
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Fig. 3. The momentum autocorrelation functions as a function of the collision 
number fit for a square system consisting of 100 disks, (a) The initial damping 
behavior of the momentum autocorrelation functions, (b) The time-oscillating re- 
gion of the momentum autocorrelation functions. C x (circles) is the autocorrelation 
function of the x-component and C y (crosses) is the autocorrelation function of the 
y-component. These two components should be equal for a square system. The sys- 
tem is the same as that in Fig. 2. The line in Fig. (b) is a fit of C x to an exponentially 
damped sinusoidal function. 



mdR + {L y -2R)/3 < < R + 2(L y -2R)/3 (the grey region in Fig. 2(a)) at 
the time of the collisions. This region is chosen because the L-mode has a node 
at this end of the system, so that the longest L-mode for 8y^ N ~^ should have 
nodes at y — R and L y — 2R and an anti-node at y — L y /2, so the Lyapunov 
mode Sy^ ^ should have largest amplitude. Further, we take a local-time 

average of the spatially averaged Lyapunov vector component 5yj 2iV_1 ^ over N 
successive disk-disk collisions. This kind of local spatial-temporal average is re- 
quired to clearly observe the time-oscillation of Lyapunov modes. A square sys- 
tem of N = 100 particles is not large enough to observe clear mode structures 
in the Lyapunov vectors, as there are only about = 10 particles in each 
coordinate direction. Fig. 2 is a graph of the local spatial-temporal average 
(5y( 2N ~^) t of dyf N ^ as a function of the collision number n t . In this figure 
we can observe the time-oscillation of the longitudinal Lyapunov mode, and 
the fit of the graph to a sinusoidal function (5y ( - 2N ^) t = asm(2irn t /T Lya + (3) 
with the fitting parameters a, (3 and Ti ya . This gives its time-oscillating period 
Tiya as Ti ya 2250. Note that Ti ya is in units of collision numbers n t , and 
the real time interval of the time-oscillating period of this Lyapunov mode is 
given by rT Lya where r ^ 0.0279 is the mean free time. 

Next, we consider the time-oscillating behavior of the momentum autocor- 
relation function. For this purpose we introduce the autocorrelation func- 
tions C v (t) of the components of momentum (where 77 = x or y) and use 
the normalized expression C v (t) = C , jy (t)/(7, ) (0), in which C v (t) is denned 
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Fig. 4. The momentum autocorrelation function C x as a function of the collision 
number rit in a rectangular system consisting of 780 disks, (a) The initial damping 
of C x for < n t < 1400. Inset: The negative region of C x for 1500 < n t < 6000. (b) 
The time-oscillating part of C x for 10000 < n t < 120000. The line is a fit of C x to 
an exponentially decaying sinusoidal function. 



by C v (t) = \im T ^ +oc ^ Cjli J T ds p vj (s + t)p vj (s). Here p vj (s) is the rj- 
component of the momentum of the j'-th disk. (In numerical calculations, 
the time-average for the time s in the above definition of autocorrelation 
function is replaced by the arithmetic average of p V j(s + t)p, n j(s) after each 
collision.) Figure 3 shows the autocorrelation functions C x and C y for the 
x- and y-components of momenta, respectively, as functions of the collision 
number n t in the same system as that used in Fig. 2. The autocorrelation 
functions C x and C y coincide as they should for a square system. Fig. 3(a) 
shows the initial damping behavior of the autocorrelation functions, and Fig. 
3(b) is the negative region in which the time-oscillating behavior of the au- 
tocorrelation functions appears. We can extract the time-oscillating period 
Ttcf by fitting this region of the autocorrelation function to the function 
C x = aexp{— f3n t } + 7exp{— en t } sm(2nn t /T tc f + () (the line in Fig. 3(b)) 
with the fitting parameters a, (3, 7, e and (. This fit gives the value of the pe- 
riod for the momentum autocorrelation function as T tc f pa 1190. Comparison 
of the value T tc f with the longest time-oscillating period T^ ya for the Lyapunov 
modes gives further evidence to support the relation (1) in two-dimensional 
systems. 



3 A Rectangular System with 780 Disks 

A system of 100 particles is not sufficient to investigate the Lyapunov modes in 
two-dimensions as only a few steps of Lyapunov spectrum are observed, even 
for hard-wall boundary conditions. A larger two-dimensional system of 780 
particles with hard-wall boundary conditions has already been investigated in 
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Ref. [9] and the Lyapunov modes calculated. The system parameters used 
were: particle mass M — 1, disk radius R = 1/2, aspect ratio L y /L x = 
0.867, and particle number density N{2R) 2 / (L x L y ) = 0.8. In this system 
the observed time-oscillating period for the Lyapunov mode is T Lya r pa 12.80, 
namely T Lya pa 74000 collisions numbers with the mean free time of r pa 
0.000173. 

Figure 4 is the normalized autocorrelation function C x of the x-component 
of the momentum as a function of the collision number n t . Here, the system 
shape is rectangular and L x > L y . This means that the time-oscillating period 
of C x is larger than that for C y , and should correspond to the T^ ya observed 
in Ref. [9]. The autocorrelation function C x , normalized by its initial value, 
shows an initial exponential damping, as shown in Fig. 4(a), and after that 
a negative region appears, as in Fig. 3 for the N = 100 case. The time- 
oscillating behavior of C x appears after the negative region of C x , as shown 
in Fig. 4(b). To extract the time-oscillating period T tc f for C x we fitted this 
time-oscillating region of the autocorrelation function to the function C x = 
aexp{— / 9n t }+7exp{— en t } sm(27in t /T tc f+C) with the fitting parameters a, (3, 
7, e and (. This fit gives us a value of T tc f ~ 40800 for time-oscillating period of 
C x . This result again confirms that the relation (1) is satisfied approximately 
for two-dimensional systems. 

We emphasize that this time-oscillating behavior of the autocorrelation func- 
tion C x appears much later than its initial damping, and with a very small 
amplitude. (The above fit gives 7 ~ 0.00127 as the amplitude of the time- 
oscillation of C x shown in Fig. 4(b), which is much smaller than the initial 
value 1 for C x .) Besides, the time-oscillation of the autocorrelation function 
has the time-scale much longer than its decay time. This means that very long 
numerical calculations are required to obtain the time-oscillating behavior of 
C x in such a large system, and we calculated more than 5 x 10 7 collisions to 
get the data presented in Figs. 3(b) and 4(b). These are much longer than the 
simulations required to calculate Lyapunov spectra. 



4 Conclusion 

In conclusion, this paper gives further numerical evidence that the relation 
T Lya pa 2T tc f for the smallest time-oscillating period T Lya of the Lyapunov 
mode and the time-oscillating period T tc f of the momentum autocorrelation 
function is satisfied in fully two-dimensional hard-disk systems. This supports 
our previous results reported in Ref. [7,8], adding numerical evidence that the 
relation is independent of the spatial dimension of the system. 

As a remark, in this paper we used a simple fit of time-oscillation of the 
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autocorrelation function and Lyapunov modes to functions including a sinu- 
soidal function of time in order to extract their time-oscillating periods, as was 
done for the quasi-one-dimensional systems [7,8]. However, the amplitude of 
the time-oscillating of momentum autocorrelation function is smaller for two- 
dimensional systems than for quasi-one-dimensional systems. Considering the 
Fourier transformation of momentum autocorrelation function and Lyapunov 
modes may be more efficient technique to extracting the periods. 

The authors of this paper appreciate the financial support of the Japan Society 
for the Promotion of Science. 
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